Contexte

Ce tutoriel a été développé dans le cadre du 28e colloque du Chapitre Saint-Laurent, présenté le 29 Mai 2024 à l’Université Laval, Québec, Canada.

Nous verrons les concepts de base qui vous permettront d’implémenter des modèles hiérarchiques avec le logiciel R, et d’interpréter les résultats en utilisant l’inférence Bayésienne.

Les objectifs d’apprentissage sont:

1. De se familiariser avec la théorie sur les modèles linéaires hiérarchiques
2. De reconnaître les variantes des modèles linéaires ((H)(G)LM)
3. D’implémenter un modèle hiéarchique avec R
4. D’interpréter les résultats avec l’inférence Bayésienne

librairies

Voici la liste des librairies que nous allons utiliser pour l’atelier. Assurez-vous de les avoir téléchargées.

# Pour manipuler les données
library(data.table)

# Pour estimer les paramètres des modèles
library(brms)

# Pour diagnostiquer les modèles bayésiens
library(bayesplot)

# Pour nos graphiques
library(ggplot2)
library(ggpubr)
library(viridis)

# Pour générer de jolies tables
library(knitr)

Jeu de données

Le jeu de données que nous utiliserons pour l’atelier se nomme ADORE et provient de l’article suivant:

Schür, C., Gasser, L., Perez-Cruz, F., Schirmer, K., Baity-Jesi, M. (2023). A benchmark dataset for machine learning in ecotoxicology. Scientific Data, 10:718. https://doi.org/10.1038/s41597-023-02612-2

Vous pouvez trouver les données originales sur le portail du Eawag Research Data Institutional Collection.

Le jeu de données ADORE contient une myriade d’informations sur des expériences écotoxicologiques effectuées sur des espèces de poisson, de crustacés, et d’algues. L’objectif des données ADORE est de concerter des informations sur la toxicité aquatique accrue de différents composés chimiques sur des espèces pour faire de la modélisation prédictive.

Le théorème de Bayes

En inférence Bayésienne, on se sert de nos connaissances à priori ainsi que des données pour faire de l’inférence. Nous pouvons représenter ceci avec le théorème de Bayes qui va comme suit:

\[P(H \mid E) = \frac{P(H) \ P(E \mid H)}{P(E)}\]

Pour résumer l’équation:

• \(P(H)\) est la probabilité d’observer notre hypothèse, c’est la distribution à priori
• \(P(E)\) est la probabilité d’observer de nouvelles preuves sous toutes les hypothèses possibles
• \(P(H \mid E)\) est la probabilité d’observer notre hypothèse selon l’évidence, c’est la distribution à postériori
• \(P(E \mid H)\) est la probabilité d’observer l’évidence si notre hypothèse est vraie, c’est la distribution de vraisemblance

Plus intuitivement, on peut aussi l’illustrer comme suit:

Illustration du théorème de Bayes issue de Yanagisawa et al. 2019

Ainsi, l’inférence Bayésienne est une approche où l’on adapte nos connaissances (à posteriori) à partir de nos hypothèses (connaissances à priori) et les preuves que nous avons.

La distribution à posteriori est donc l’élément central que nous cherchons à estimer. Rappelez-vous que chaque fois que vous voyez le terme “à posteriori” dans ce tutoriel, nous faisons référence à la distribution d’un paramètre d’intérêt estimé par le modèle.

En prime, voici un lien vers une application qui vous permet de jouer avec ces trois distributions pour mieux comprendre comment les connaissances à priori et les preuves modifient la distribution à posteriori:

https://micl.shinyapps.io/prior2post/

Méthode d’estimation

Il existe de nombreuses façons d’approximer la distribution à posteriori. Dans ce tutoriel, nous utiliserons la méthode de Monte Carlo Hamiltonienne (HMC) et l’échantillonneur No-U-Turn (NUTS) pour estimer la distribution à posteriori. Ce sont des algorithmes spécifiques dans la grande famille d’estimation de Monte Carlo par chaîne de Markov (MCMC).

Pour ce faire nous effectuerons tous nos modèles avec la librairie brms, une interface R pour le langage probabiliste STAN. La librairie brms utilise l’algorithme NUTS pour estimer les paramètres. Je vous recommande vivement de lire les vignettes de la librairie pour mieux comprendre les capacités et le fonctionnement de brms.

Les avantages

L’inférence Bayésienne comporte certains avantages par rapport à l’approche statistique fréquentiste.

• Elle met l’accent sur la distribution et l’incertitude au lieu d’estimés uniques
• Elle permet une interprétation plus naturelle (probabiliste)
• On peut facilement critiquer les modèles

Pour aller plus loin

La théorie sous-jacente au théorème de Bayes dépasse largement le cadre de ce tutoriel, donc veuillez consulter les ressources pertinentes dans la section Références. Une très bonne introduction que je ne peux recommander assez est l’excellent livre “Statistical Rethinking” de Richard McElreath. Richard McElreath publie également l’ensemble du contenu du livre sous forme de vidéos sur YouTube.

Description générale

Pour l’atelier, nous utiliserons une version des données nettoyées et prêtes pour l’analyse qui se trouvent sur ce dépôt GitHub dans le dossier data. Ce dépôt a été préparé spécifiquement pour l’atelier.

Voici les spécifications:

• Le jeu de données servira à faire des modèles linéaires généralisés (GLM) et des modèles linéaires généralisés hiérarchiques (HGLM) univariés
• Nous appliquerons les modèles aux données d’animaux seulement. Nous appelerons le jeu de données anim

Importons les données à partir du dépôt GitHub de l’atelier:

github <- "https://raw.githubusercontent.com/quantitative-ecologist"
depot <- "Atelier-CSL-EcotoQ"
dossier <- "main/data"

anim <- fread(
  file.path(github, depot, dossier, "animal_dat.csv"),
  stringsAsFactors = TRUE
)

Variables

Une description détaillée des variables du jeu de données est présentée en Annexe de l’article de Schür et al. (2023).

Dans le cadre de cet atelier, nous utiliserons un échantillon de ces variables.

Distribution de la concentration effective

Commençons par évaluer la distribution de la variable réponse, qui est le log de la concentration effective moyenne. Nous allons faire nos graphiques avec la librairie ggplot2.

ggplot(data = anim) +
  geom_histogram(
    aes(x = result_conc1_mean_log),
    color = "black", fill = "#440154"
  ) +
  ylab("Compte") +
  xlab("log(concentration moyenne (mg/L))") +
  theme_classic(base_size = 14) +
  theme(panel.grid = element_blank())

Concentration effective selon l’espèce

Vérifions si les concentrations effectives varient selon l’espèce.

ggplot(data = anim) +
  geom_boxplot(
    aes(
      x = result_conc1_mean_log,
      y = tax_gs,
      fill = tax_group
    ),
    color = "black"
  ) +
  ylab("Espèce") +
  xlab("log(concentration moyenne (mg/L))") +
  scale_fill_manual(values = c("#21918c", "#440154")) +
  scale_x_continuous(
    breaks = seq(-8, 8, 4),
    limits = c(-9, 9)
  ) +
  theme_classic(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "top"
  )

Concentration effective selon l’exposition

Finalement, on peut aussi se douter que la concentration effective pourrait changer selon le type d’exposition.

ggplot(data = anim) +
  geom_boxplot(
    aes(
      x = result_conc1_mean_log,
      y = test_exposure_type,
      fill = tax_group
    ),
    color = "black"
  ) +
  ylab("Type d'exposition") +
  xlab("log(concentration moyenne (mg/L))") +
  scale_fill_manual(values = c("#21918c", "#440154")) +
  scale_x_continuous(
    breaks = seq(-8, 8, 4),
    limits = c(-9, 9)
  ) +
  theme_classic(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    legend.position = "top"
  )

Sommaire de l’exploration de données

Les vérifications nous indiquent qu’il semble y avoir une structure hiérarchique dans nos données. En termes simples, cela signifie qu’il y a de la variation dans les concentrations effectives entre les différents groupes que nous avons vu, soit, les espèces ainsi que le type d’exposition.

Les modèles hiérarchiques nous permettront de modéliser cette variation afin d’en tenir compte.

Les bases

L’objectif d’un modèle linéaire est de prédire une variable \(Y\) tout en minimisant l’écart entre les observations et les prédictions (\(\hat{Y}\)).

Au travers votre parcours, vous avez sûrement vu différents types de modèles où on utilise des variables explicatives \(X_i\) pour prédire une variable \(Y\) (ex. ANOVA, régression, ANCOVA).

Par exemple, une régression linéaire peut s’écrire avec cette équation:

\[Y = X\beta+\epsilon\]

où:

• \(Y\) est le vecteur des valeurs observées (c.-à.-d. la variable réponse)
• \(X\) est la matrice de variables explicatives (incluant l’ordonnée à l’origine)
• \(\beta\) est le vecteur de coefficients
• \(\epsilon\) est le vecteur des résidus

Visuellement, on peut représenter ce modèle avec la figure suivante:

Pour interpréter les résultats d’un modèle linéaire, il est essentiel que celui-ci respecte ces conditions d’application:

1. Linéarité: la relation entre \(Y\) et \(X\) est linéaire.
2. Indépendance: les observations sont indépendantes l’une de l’autre.
3. Homoscédasticité: la variance des résidus est homogène le long des valeurs prédites.
4. Normalité des résidus: les résidus suivent une distribution Gaussienne.

Limites des modèles linéaires

Une des limites majeures des modèles linéaires classiques est qu’ils ne peuvent s’ajuster à des variables réponse pour lesquelles les résidus suivent une distribution non-Gaussienne.

Cette situation est particulièrement commune lorsqu’on analyse des données biologiques, qui ont tendance à suivre d’autres types de distribution telle que de type Poisson ou binomiale.

Dans la section suivante, nous verrons comment les modèles linéaires généralisés nous permettent de palier à ce problème en modélisant spécifiquement la nature des résidus d’une variable \(Y\).

Familles de distribution

Il existe une multitude de familles de distribution qui décrivent les variables qu’on peut mesurer. Ici, l’important est de retenir que chaque distribution a des paramètres propres qui décrivent sa forme. Les GLM utilisent ces paramètres pour l’estimation des coefficients, et nous les utiliserons pour l’interprétation de nos modèles.

Voici trois exemples des familles les plus communes.

Fonction de lien

Nous avons vu trois exemples de distributions souvent rencontrées avec des données biologiques.

Dans un GLM, afin que celui-ci puisse estimer la relation entre \(Y\) et \(X\), on doit donc:

1. Spécifier la distribution des résidus
2. Spécifier la fonction de lien pour les valeurs prédites

Cette deuxième étape consiste à linéariser la relation entre \(Y\) et \(X\). Il s’agit donc de transformer les valeurs prédites \(\hat{Y}\) pour avoir une relation linéaire.

Mathématiquement, le prédicteur linéaire dans un GLM est défini comme:

\[\eta = X\beta\]

où:

• \(\beta\) sont les coefficients
• \(X\) est la matrice des variables

La fonction de lien transforme le prédicteur linéaire à l’échelle appropriée pour \(Y\) tel que:

\[g(\mu) = \eta\]

où:

• \(\mu\) est la moyenne prédite de \(Y\) en fonction des prédicteurs \(X\)
• \(g(.)\) est la fonction de lien

Voici une table qui résume les informations relatives à la fonction de lien pour les familles que nous avons vu.

Préparation des données

Commençons par évaluer la structure des variables avant de faire le modèle.

str(anim[
  , .(
    result_conc1_mean_binary, test_exposure_type,
    tax_lh_amd, media_temperature_mean
  )
])

## Classes 'data.table' and 'data.frame':   28999 obs. of  4 variables:
##  $ result_conc1_mean_binary: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ test_exposure_type      : Factor w/ 4 levels "F","NR","R","S": 3 3 3 3 3 3 3 3 3 3 ...
##  $ tax_lh_amd              : num  5475 5475 5475 5475 5475 ...
##  $ media_temperature_mean  : num  19 19 19 19 19 19 NA NA NA NA ...
##  - attr(*, ".internal.selfref")=<externalptr>

Les données sont dans le format adéquat pour la modélisation.

Avans de faire le modèle, nous allons normaliser les variables explicatives continues en les ramenant à la même échelle. Ceci peut faciliter les calculs des coefficients et simplifier l’interprétation des paramètres.

# Fonction de normalisation en score Z.
standardize <- function(x) {
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

# Créer deux nouvelles colones
anim[
  , c("Ztax_lh_amd", "Zmedia_temperature_mean") := lapply(
    .SD,
    standardize
  )
  , .SDcols = c("tax_lh_amd", "media_temperature_mean")
]

Construction du modèle

Dans cet exemple, on s’intéresse à estimer comment la toxicité de la concentration effective varie selon le type d’exposition, la longévité moyenne de l’espèce, et la température du médium durant l’expérience.

Voici l’équation de notre modèle:

\[\eta = X\beta + Z\gamma + 1\beta_0 + \epsilon\]

où:

• \(\eta\) est le prédicteur linéaire
• \(X\) est la matrice de conception du type d’exposition
• \(\beta\) est le vecteur des coefficients du type d’exposition
• \(Z\) est la matrice de conception des prédicteurs linéaires (longévité et température)
• \(\gamma\) est le vecteur des coefficients du type d’exposition des prédicteurs linéaires
• \(1\) est le vecteur de 1 pour l’ordonnée à l’origine
• \(\beta_0\) est l’ordonnée à l’origine pour le niveau de référence

La fonction de lien du modèle est:

\[logit(p) = log(\frac{p}{1-p})\]

où:

\(p\) est la probabilité que la concentration effective soit toxique (1) pour chaque observation.

On estime \(p\) en utilisant l’inverse de la fonction logit:

\[p = \text{logit}^{-1}(\eta) = \frac{1}{1 + e^{-\eta}}\]

Spécifions ce modèle dans brms avec la fonction brmsformula():

form1 <- brmsformula(
  result_conc1_mean_binary ~
  test_exposure_type +
  Ztax_lh_amd +
  Zmedia_temperature_mean
)

Informations à priori

Comme nous l’avons vu dans la section sur l’inférence Bayésienne, la distribution postérieure est estimée en multipliant la distribution à priori et la distribution de vraisemblance des données (les preuves).

Par défaut, brms emploie des distributions à priori très larges qui ne réflètent pas nos hypothèses/connaissances sur notre système. Définir les connaissances à priori est donc très important en inférence Bayésienne, et je vous suggère fortement de faire ce processus. C’est ce qui différencie l’inférence Bayésienne de l’approche fréquentiste.

Ici, puisque je ne connais pas assez bien l’écotoxicologie, je spécifie des “priors” relativement informatifs suivant une distribution Gaussienne avec une moyenne de 0 et un écart-type de 1.

Cela signifie que je crois que les valeurs des paramètres devraient se situer autour de 0 avec des valeurs positives allant de +1 à -1 écart-types.

Les informations à priori à spécifier sont:

• sur l’ordonnée à l’origine globale (la moyenne du type d’exposition F)
• sur les différences de moyenne par rapport à F
• sur les pentes pour la température et la longévité

Pour vous aider à comprendre, imaginez que vous savez avec la littérature que la dose effective devrait être plus toxique avec une augmentation de la température. Ainsi, pour la pente de ce paramètre, vous pourriez spécifier une distribution Gaussienne avec une moyenne de 1 et un écart-type de 0.5. Vous diriez donc à votre modèle que vous assumez que cette valeur devrait être positive, se situant autour de 1, avec peu de variation autour de cette valeur.

priors1 <- c(
  # Informations à priori sur l'ordonnée à l'origine
  set_prior(
    "normal(0, 1)",
    class = "Intercept"
  ),
  # Informations à priori sur les différences de moyenne
  set_prior(
    "normal(0, 1)",
    class = "b",
    coef = c(
      "test_exposure_typeNR",
      "test_exposure_typeR",
      "test_exposure_typeS"
    )
  ),
  # Informations à priori sur les pentes
  set_prior(
    "normal(0, 1)",
    class = "b",
    coef = c(
      "Zmedia_temperature_mean",
      "Ztax_lh_amd"
    )
  )
)

Estimation des paramètres

parallel::detectCores()

## [1] 16

glm1 <- brm(

  # La formule spécifiée
  formula = form1,

  # La famille de distribution
  family = bernoulli(link = "logit"),

  # Les paramètres d'itérations pour le MCMC
  iter = 2000,
  warmup = 1000,
  thin = 1,
  chains = 4,

  # Paramètres pour accélérer les calculs
  # Utiliser seulement si vous avez >4 coeurs
  threads = threading(3),

  # Si vous avez 4 coeurs ou moins
  # cores = parallel::detectCores()-1,

  # Estimateur
  backend = "cmdstanr",

  # Reproductibilité
  seed = 123,

  # Paramètres de contrôle
  control = list(adapt_delta = 0.95),

  # Distributions à priori
  prior = priors1,

  # Échantillonage des distributions à priori
  # Ceci peut être utile à des fins de comparaison
  sample_prior = TRUE,

  # Données
  data = anim
)

chemin <- file.path(getwd(), "outputs")
dir.create(chemin)

saveRDS(
  object = glm1,
  file = file.path(chemin, "glm1.rds")
)

Vérifications du modèle

Maintenant que notre modèle a terminé son exécution, il est essentiel de faire des inspections pour s’assurer que les paramètres estimés sont fiables. Ces étapes d’inspections sont cruciales pour s’assurer que nos interprétations biologiques font du sens. Autrement, nous pourrions avoir de mauvaises surprises.

Commençons par extraire les échantillons des distributions à posteriori dans un data.frame pour les manipuler.

params <- as_draws_df(
  x = glm1, add_chain = TRUE,
  regex = TRUE
)

bayesplot::mcmc_trace(
  params,
  regex_pars = "b",
  np = nuts_params(glm1)
)

glm1_t <- data.table(
  round(summary(glm1, prob = 0.89, robust = TRUE)$fixed, digits = 3),
  keep.rownames = TRUE
)
setnames(glm1_t, old = "rn", new = "Parameter")

glm1_t[, c(1, 6:8)]

pp_check(glm1, ndraws = 100)

Interprétations

glm1_t <- data.table(
  round(summary(glm1, prob = 0.89, robust = TRUE)$fixed, digits = 3),
  keep.rownames = TRUE
)
setnames(glm1_t, old = "rn", new = "Parameter")

knitr::kable(glm1_t[, c(1, 2, 4, 5)])

tab1 <- conditional_effects(
  x = glm1, method = "posterior_epred",
  effects = "Zmedia_temperature_mean",
  spaghetti = TRUE, ndraws = 100
)
spag1 <- attr(tab1$Zmedia_temperature_mean, "spaghetti")

tab2 <- conditional_effects(
  x = glm1, method = "posterior_epred",
  effects = "Ztax_lh_amd",
  spaghetti = TRUE, ndraws = 100
)
spag2 <- attr(tab2$Ztax_lh_amd, "spaghetti")

params

# Transformer en data.table pour manipuler
paramsdt <- data.table(params)

# On calcule les vraies moyennes
paramsdt[, c(2:4) := paramsdt[, c(2:4)] + b_Intercept]

# Changer en format long
preds3 <- melt(
  data = paramsdt[, c(1:4)],
  variable.name = "exposure_type",
  measure = colnames(paramsdt[, c(1:4)])
)

# Ramener à l'échelle de probabilité
preds3[, value := plogis(value)]

# Calculer les médianes
medians <- preds3[, .(median = median(value)), by = exposure_type]

ggplot(
  data = preds3,
  aes(x = value, fill = exposure_type)
) +
  geom_density(alpha = 0.5, colour = "black") +
  geom_vline(
    data = medians,
    aes(xintercept = median, color = exposure_type),
    linetype = "dashed"
  ) +
  scale_fill_manual(
    labels = c(
      "Flux", "NS",
      "Renouvellement", "Statique"
    ),
    values = c("#440154", "#3d4d8a", "#24878e", "#40bd72")
  ) +
  scale_colour_manual(
    labels = c(
      "Flux", "NS",
      "Renouvellement", "Statique"
    ),
    values = c("#440154", "#3d4d8a", "#24878e", "#40bd72")
  ) +
  #scale_x_continuous(breaks = seq(0.25, 1, 0.25), limits = c(0.25, 1)) +
  xlab("\nProbabilité d'être plus toxique") +
  ylab("Densité\n") +
  labs(fill = "Exposition", colour = "Exposition") +
  theme_classic(base_size = 14) +
  theme(legend.position = "top")

preds1 <- ggplot() +
  geom_point(
    data = anim,
    aes(
      y = result_conc1_mean_binary,
      x = Zmedia_temperature_mean
    ),
    alpha = 0.05
  ) +
  geom_line(
    data = spag1,
    aes(
      x = Zmedia_temperature_mean,
      y = estimate__,
      group = sample__
    ),
    linewidth = 1, alpha = 0.05,
    colour = "#440154"
  ) +
  ylab("Toxicité\n") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  xlab("\nTempérature moyenne normalisée") +
  theme_classic(base_size = 14)

preds2 <- ggplot() +
  geom_point(
    data = anim,
    aes(
      y = result_conc1_mean_binary,
      x = Ztax_lh_amd
    ),
    alpha = 0.05
  ) +
  geom_line(
    data = spag2,
    aes(
      x = Ztax_lh_amd,
      y = estimate__,
      group = sample__
    ),
    linewidth = 1, alpha = 0.05,
    colour = "#440154"
  ) +
  ylab("Toxicité\n") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  xlab("\nLongévité moyenne") +
  theme_classic(base_size = 14)

ggarrange(preds1, preds2)

Complexité des données biologiques

En explorant nos données, nous avons vu qu’il semble y avoir une structure hiérarchique, où la concentration effective varie entre les espèces.

Ceci est un exemple concret démontrant comment les données biologique sont souvent complexes, et structurées à différentes échelles.

Reprenons l’exemple de la concentration effective en fonction de la température du médium. Cette fois, on ajoute les espèces dans la figure. Ici, chaque couleur correspond à un taxon différent.

À partir de ce graphique, on peut se poser deux questions:

1. Est-ce que les taxons ont la même concentration effective moyenne?
2. Est-ce que la relation entre la concentration effective et la température change selon les taxons?

Pour la première question, notre réflexe pourrait être de faire le même GLM que nous avons fait à la section précédente, et d’ajouter le taxon comme variable catégorique.

Pour la seconde question, on pourrait spécifier un modèle de type ANCOVA, en modélisant une interaction entre la température et le taxon.

Pour répondre aux deux questions, on aurait ce modèle:

result_conc1_mean_log ~
  test_exposure_type +
  Ztax_lh_amd +
  Zmedia_temperature_mean +
  tax_gs +
  Zmedia_temperature_mean:tax_gs

Toutefois, on se retrouve avec un gros problème ici. Visualisons-le:

Nous avons 105 espèces, et donc 105 droites à estimer. Cela signifie que le modèle doit estimer:

• 105 ordonnées à l’origine
• 105 pentes

C’est donc 210 paramètres de + à estimer! Si on utilisait l’approche fréquentiste pour l’analyse de données, ça nous coûterait cher en degrés de liberté. Ceci est moins problématique pour l’approche Bayésienne car on n’utilise pas la valeur p.

Fréquentiste ou Bayésien, le plus important, c’est qu’en estimant 210 paramètres, la sortie de la fonction summary() serait cauchemardesque. Pensez-y! Plus haut, nous avions un Intercept_b qui utilisait le type d’exposition F comme référence. Ensuite, le modèle calculait la différence entre les autres types d’exposition et F. Je vous laisse imaginer celà avec 105 taxons…

L’avantage des modèles hiérarchiques

Les modèles hiérarchiques nous permettent de palier aux problèmes que nous avons vu plus haut:

• permettent de séparer des effets populationnels d’effets groupés
• permettent de modéliser la structure hiérarchique de nos données en réduisant le nombre de paramètres à estimer
• exploitent les paramètres populationnels pour estimer les paramètres de groupe même quand la taille d’échantillon est faible pour un groupe

De plus, rappelez vous qu’une des conditions de base dans les modèles linéaires est que les échantillons doivent être indépendants. Dans l’exemple que nous avons vu, nous avons plusieurs observations par espèce, ce qui contrevient à cette condition d’indépendance.

Les modèles hiérarchiques permettent donc de modéliser des facteurs associés à la structure de vos expériences. Voici quelques exemples:

• On peut intégrer les quadrats, ou les sites
• On peut ajouter les individus si des données sont prises à répétition
• Permettent de considérer l’aspect temporel dans les données

Nous verrons comment ces principes fonctionnent dans la prochaine section avec les concepts d’effets fixes et d’effets aléatoires.

Effets fixes et aléatoires

Reprenons l’exemple du modèle que nous avons construit précédemment.

On s’intéresse toujours à estimer comment la concentration effective varie selon le type d’exposition, la longévité moyenne de l’espèce, et la température du médium durant l’expérience. On ajoute toutefois les espèces comme effet aléatoire.

Ici, on veut donc savoir si les espèces varient dans leur concentration effective moyenne, et si la relation entre celle-ci et la température change d’une espèce à l’autre.

Voici l’équation de notre modèle:

\[\eta = X\beta + Z\gamma + U_1\delta_1 + U_2\delta_2 + 1\beta_0 + \epsilon\]

où les paramètres ajoutés sont:

• \(U_1\) est la matrice de conception pour les ordonnées à l’origine aléatoires
• \(\delta_1\) est le vecteur des ordonnées à l’origine aléatoires des espèces
• \(U_2\) est la matrice de conception pour les pentes aléatoires
• \(\delta_2\) est le vecteur des pentes aléatoires des espèces

Les effets fixes, représentés par les coefficients \(\beta\) et \(\gamma\), modélisent les effets moyens des variables explicatives sur la variable réponse. Ils sont constants pour toutes les unités observées (c.-à.-d., les espèces par exemple). Ce sont les effets que nous avons estimés précédemment dans notre GLM. Celà ne change donc pas. Les effets fixes peuvent être des variables catégoriques ou des variables continues

Les effets aléatoires, quant à eux, représentés par les coefficients \(\delta_1\) et \(\delta_2\), modélisent la variation entre les différents niveaux d’une variable catégorique. Dans le contexte de notre exemple, on intègre les différences entre les espèces en termes de sensibilité à la concentration en ajoutant un effet aléatoire. Les effets aléatoires sont strictement des variables catégoriques

En ajoutant cet effet aléatoire, le modèle peut capturer une part supplémentaire de la variation de \(Y\) qui ne peut pas être expliquée par les effets fixes. Cela permet de mieux ajuster le modèle aux données et d’améliorer les prédictions.

Les types de modèles hiérarchiques

Il existe deux types de modèles hiérarchiques correspondant aux effets aléatoires estimés.

1. Pour le premier type, on assume que la relation entre \(Y\) et \(X\) est la même pour l’ensemble de nos groupes, mais que la moyenne de \(Y\) varie entre les groupes. C’est donc un modèle où l’on estime une moyenne pour chacun des niveaux de notre facteur d’intérêt, mais pour lequel la relation entre \(Y\) et \(X\) est la même que l’effet fixe pour l’ensemble des niveaux du facteur. Ce type de modèle hiérarchique se nomme modèle d’ordonnées à l’origine aléatoires.

2. Le second type de modèle est construit à partir de la mécanique du modèle d’ordonnées à l’origine aléatoires. C’est le modèle de pentes aléatoires. Comme le nom l’indique, c’est un modèle où chaque niveau du facteur peut maintenant avoir sa propre relation entre \(Y\) et \(X\).

Comme nous l’avons soulevé plus haut, un avantage majeur des modèles hiérarchiques est qu’ils réduisent le nombre de paramètres à estimer. Au début de cette section, nous avons vu que si nous faisions une ANCOVA, nous aurions 210 paramètres à estimer. En exploitant les effets aléatoires nous réduisons ce nombre de paramètres à 2 seulement!

La raison est qu’au lieu d’estimer la différence d’ordonnée à l’origine pour chaque niveau du facteur ainsi que la différence de pente, un modèle hiérarchique va plutôt calculer chaque valeur pour chaque niveau (c.-à-d., moyenne et pente), et ensuite estimer une variance.

On a donc:

• la variance entre les ordonnées à l’origine (c.-à-d., les moyennes des niveaux du facteur)
• la variance entre les pentes (c.-à-d., la relation entre \(Y\) et \(X\) pour chaque niveau du facteur)

Ce que la variance nous dit:

• une variance faible entre les ordonnées à l’origine = la moyenne de \(Y\) de chaque espèce est semblable
• une variance faible entre les pentes = la relation entre \(Y\) et \(X\) pour chaque espèce est semblable

On assume que les ordonnées à l’origine et les pentes suivent une distribution normale telle que:

\[ \begin{bmatrix} \delta_1 \\ \delta_2 \end{bmatrix} \sim \text{MVN} \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \Omega \right) \]

où \(\Omega\) est une matrice de covariance telle que:

\[ \Omega = \begin{bmatrix} \text{var}(\delta_1) & \text{cov}(\delta_1, \delta_2) \\ \text{cov}(\delta_2, \delta_1) & \text{var}(\delta_2) \end{bmatrix} \]

Inspections additionnelles

While the model verifications that we just did are probably among the most important ones, they cover only a small portion of the checks that you should perform to ensure that your model has been correctly specified. To go further, here is a list of useful links for additional checks:

• MCMC diagnostics with the bayesplot package
• Additional posterior predictive checks with the bayesplot package
• Additional residuals inspections with the tidybayes package. For instance, you could inspect the relationship between the residuals and your fixed effects to see if the model failed at describing some features. This will be important for biological interpretations.
• Leave-one-out cross-validation using Pareto-smoothed importance sampling for predictive performance and model comparisons. For important references consult the loo FAQ, the roaches model example, and the loo package vignette

Moreover, note that MCMC estimation can come with different convergence problems. Those problems can arise from different processes, but brms (and other packages using MCMC) should print warnings outlining the potential problems after the model has finished running. Here is a useful link describing each of those problems and their solutions.